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The decoherence of a quantum system S coupled to a quantum environment E is considered. For states 
chosen uniformly at random from the unit hypersphere in the Hilbert space of the closed system S + E we derive 
a scaling relationship for the sum of the off-diagonal elements of the reduced density matrix of 5 as a function 
of the size De of the Hilbert space of E. This sum decreases as 1/i/Dg as long as De 2> l. This scaling 
prediction is tested by performing large-scale simulations which solve the time-dependent Schrodinger equation 
for a ring of spin-l/2 particles, four of them belonging to S and the others to E. Provided that the time evolution 
drives the whole system from the initial state toward a state which has similar properties as states belonging 
to the class of quantum states for which we derived the scaling relationship, the scaling prediction holds. For 
systems which do not exhibit this feature, it is shown that increasing the complexity (in terms of connections) 
of the environment or introducing a small amount of randomness in the interactions in the environment suffices 
to observe the predicted scaling behavior. 

PACS numbers: 03.65.Yz, 75.lO.Jm, 75.l0.Nr, 05.45.Pq 



I. INTRODUCTION 

Decoherence of a quantum system S interacting with a 
quantum environment E is of importance for two reasons. 
First, decoherence of S is the primary requirement for S to 
relax to a state described by a canonical ensemble at a certain 
temperature JT)]. Second, decoherence is arguably the largest 
impediment for practical, realizable quantum computers fl2|]. 

The large interest in technological areas like spintron- 
ics, quantum computing and quantum information process- 
ing have stimulated the theoretical research of quantum dy- 
namics in open and closed interacting systems. Besides this 
more application driven interest there persists the fundamental 
and still unanswered question under which conditions a finite 
quantum system reaches thermal equilibrium and how this can 



be derived from dynamical laws. 

On the one hand there exists a variety of studies explor- 
ing the microcanonical thermalization in an isolated quantum 
system On the other hand there exist various studies 

investigating the process of canonical thermalization of a sys- 
tem coupled to a (much) larger system I7l-fl3l and of two 
finite identical quantum systems prepared at different temper- 
atures (30. 

In previous work iflrj [l7ll . we numerically demonstrated 
that a quantum system interacting with an environment at high 
temperature relaxes to a state described by the canonical en- 
semble. In this paper we focus on investigating the dynamic 
properties of the decoherence of a quantum system S, being 
a subsystem of the whole system S + E. We do this both 
with a theoretical prediction and by simulating the dynamics 
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of a relatively large system S + E of spin-1 /2 particles using 
a time-dependent Schrodinger equation (TDSE) solver lfl8ll . 
In particular, we investigate the scaling of the degree of deco- 
herence of S with the size of E, keeping the size of S fixed. 
Based on similar arguments as given in Ref. Il9ll . we find that 
the degree of decoherence of S decreases as 1 / \/De, where 
De is the dimension of the Hilbert space of the environment if 
the state of the whole system is chosen uniformly at random 
from the unit hypersphere in the Hilbert space. In this paper, 
we denote states chosen uniformly at random from the unit 
hypersphere in the Hilbert space of the whole system by "X" 
and of the environment by "Y". 

We also address the question under what circumstances the 
whole system evolves to a state which has the same degree of 
decoherence as a state "X". In particular we study the case 
in which the initial state of S + E is a direct product of the 
state Itltl) °f S and a state "7" of E. If the initial state of the 
whole system S + E is slightly different from a given state "X", 
the dynamics may drive the whole system into a state which 
is very different from the given state "X", but which is of a 
similar type. We investigate through our simulations when the 
dynamics plays an important role in the decoherence in that it 
can drive S+E to a state "X" by introducing small world bond 
connections in E and/or between S and E and by introducing 
randomness in the interaction strengths of the environment. 

The paper is organized as follows. In Section II our theo- 
retical results for the scaling of the decoherence of S are pre- 
sented, together with details of the one-dimensional ring of 
spin-1 /2 particles which we simulate to better understand the 
scaling prediction. Sections III-V contain results for the one- 
dimensional rings under study. In particular we look at the ef- 
fect of adding additional bonds (Small World Bonds, SWBs) 
between the system and environment spins and/or between en- 
vironment spins only (Section IV) and of randomness in the 
interaction strengths of the Hamiltonian of the environment 
(Section V). Section VI contains our conclusions and a dis- 
cussion of our results. 



II. THEORY, MODEL, AND METHODS 

The time evolution of a closed quantum system is governed 
by the time-dependent Schrodinger equation (TDSE) ll20[|2lll . 
If the initial density matrix of an isolated quantum system is 
non-diagonal then, according to the time evolution dictated 
by the TDSE, it remains non-diagonal. Therefore, in order 
to decohere the system S, it is necessary to have the system 
S interact with an environment E, also called a heat bath or 
spin bath if the environment is composed of spins. Thus, the 
Hamiltonian of the whole system S + E takes the form 

H = H S + H E + H SE , (1) 

where H$ and He are the system and environment Hamilto- 
nian respectively, and H$ E describes the interaction between 
the system and environment. In what follows, we first describe 
the general theory that leads to the scaling of the decoherence 
of the system S with the size of E and S. We then describe 



in detail the spin- 1/2 Hamiltonians we have simulated to pro- 
vide a case study for this scaling. 

A. Time evolution 

A pure state of the whole system S + E evolves in time ac- 
cording to (in units of h = 1 ) 

As d e 

m)) = e- m |¥(0)> = £ £ c(i,p,t) \i,p) , (2) 
i=i p=\ 

where the set of states {|z',/?)} denotes a complete set of or- 
thonormal states in some chosen basis, and D$ and De are the 
dimensions of the Hilbert spaces of the system and the envi- 
ronment, respectively. We assume that D$ and D E are both 
finite. 

The spin Hamiltonian H models a system with Ng spin- 
1/2 particles and an environment with Ne spin- 1/2 particles. 
Thus, Ds = 2 Ns and De = 2 Ne . The whole system S + E con- 
tains N = Ns + Ne spin-1 /2 particles and the dimension of its 
Hilbert space is D = D^De- In our simulations we use the 
spin-up - spin-down basis. Numerically, the real-time prop- 
agation by e~" H is carried out by means of the Chebyshev 
polynomial algorithm 1I22I - I25I1 . thereby solving the TDSE for 
the whole system starting from the initial state ^(O)). This 
algorithm yields results that are very accurate (close to ma- 
chine precision), independent of the time step used lfl8ll . 

B. Computational aspects 

Computer memory and CPU time severely limit the sizes 
of the quantum systems that can be simulated. The required 
CPU time is mainly determined by the number of operations 
to be performed on the spin-1 /2 particles. The CPU time does 
not put a hard limit on the simulation. However, the memory 
of the computer does severely limit which system sizes can 
be calculated. The state \ x ¥) of a A^-spin-1/2 system is rep- 
resented by a complex-valued vector of length D = 2 N . In 
view of the potentially large number of arithmetic operations, 
it is advisable to use 13-15 digit floating-point arithmetic 
(corresponding to 8 bytes for a real number). Thus, to repre- 
sent a state of the quantum system of N spin-1 /2 particles on 
a conventional digital computer, we need a least 2 N+4 bytes. 
Hence, the amount of memory that is required to simulate a 
quantum system with N spin-1 /2 particles increases exponen- 
tially with N. For example, for N = 24 (N = 36) we need at 
least 256 MB (1 TB) of memory to store a single arbitrary 
state I 1 ?) . In practice we need three vectors, memory for com- 
munication buffers, local variables and the code itself. 

The elementary operations performed by the computational 
kernel are of the form I 1 ?) <— U I^P) where U is a sparse uni- 
tary matrix with a very complicated structure (relative to the 
computational basis). Inherent to the problem at hand is that 
each operation U affects all elements of the state vector I 1 ?) in 
a nontrivial manner. This translates into a complicated scheme 
for accessing memory, which in turn requires a sophisticated 
MPI communication scheme (26ll . 
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C. Reduced density matrix 

The state of the quantum system S is described by the re- 
duced density matrix 



p(t)=Tr E p (f), 



(3) 



where p (?) is the density matrix of the whole system S + E at 
time t and Trg denotes the trace over the degrees of freedom 
of the environment. In terms of the expansion coefficients 
c(i,p,t), the matrix element (ij) of the reduced density ma- 
trix reads 

D E D E 

P,j(t) = Tr £ £ £c*(/,<7,f)c(;>,f)l;»('\<7l 

p=\q=\ 



E c*(i,p,t)c(j,p,t) . 



(4) 



We characterize the degree of decoherence of the system by 

(7(f) 



Dc-1 Dc 



\ E E Mo I 



(5) 



=1+1 



matrix p in the representation that diagonalizes Hg. Clearly, 
(7(f) is a global measure for the size of the off-diagonal terms 
of the reduced density matrix in the representation that diago- 
nalizes H$. If a(t ) = the system is in a state of full decoher- 
ence (relative to the representation that diagonalizes H$). 



D. Scaling property of a 

We can prove a scaling property of a by assuming that the 
final state of the whole system is a state "X", a state that is 
picked uniformly at random from the unit hypersphere in the 
Hilbert space. The wave function of the whole system reads, 



D s D E 

E E c„ P 

i= 1 p= 1 



-P J : 



(6) 



where j ^ } Ep E ^ ^ j^j is the set of eigenvectors of H$ 
(He), and the real and imaginary parts of Q tP are real ran- 
dom variables. The derivation of the scaling behavior fol- 
lows Refjll. In particular Eqs. (A8), (A12) and (A23) 
of Ref. lfl9ll are used. We introduce the following short- 
hand notation for the sum over the off-diagonal elements, 



K ij = Efil £/il ( l ~ S ij) K ij fOT an y K 'j< Where $7 is the 
where py(f) is the matrix element of the reduced density Kronecker delta function. The expectation value is given by 
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c>l \Cj, P \ )=LL^ Db{DsDe+1] 



D S D E + 1 De + jL 



(7) 



where E(-) denotes the expectation value with respect to the 
probability distribution of the random variables C;. p . Equa- 
tion <JVj> does not require any condition on the Hamiltonian 
Eq. (fTJ. For example, if He is composed of two or more en- 
vironments that do not couple to each other, but only interact 
with the system, in Eq. (O De is the product of the sizes of 
the Hilbert spaces of all the environments. In addition, Eq. (|7]i 
does not impose any requirement on the geometry. 

From Eq. (|7]i it follows that for any fixed value of Ds > 1 
and De ^> 1, (7 scales as 

V2 V V2V D s D E + l V2D^ 

Therefore, if the size of the system S is fixed (which is the case 
considered in this paper), a decreases as 1 j \J~De for large De- 
Hence, for a spin-1 /2 system a should decrease as 2~ N£ I 2 for 
large N E . 



For fixed Ds > 1, it follows from Eq. (0 that the environ- 
ment does not have to be verylarge for Eq. (|8) to hold, which 
is in agreement with Ref. 12711 . Nevertheless, the existence 
of an environment is crucial. If there is no environment, then 
the <7 approaches to a constant (see Appendix[A}, even if the 
whole system is initially in a state "X". 



E. Model and method 

For testing the predicted scaling of Eq. © we simulate sys- 
tems of spin- 1/2 particles. For studying the time evolution 
of the whole system S + E, we consider a general quantum 
spin-1 /2 model defined by the Hamiltonian of Eq. (Q~|i where 

N s - 1 N s 

Hs = - E E E W^l w 

<'=1 j=i+\ a=x.y,z 
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FIG. 1. (Color online) An example of a spin system used in the 
simulations. The N$ = 4 system spin-1/2 particles are colored light 
gray (cyan), and the Ne = 18 environment spin-1/2 particles are col- 
ored dark gray (red). The thin black segments show the connections 
for a one-dimensional ring, which are the only bonds (interactions) 
present in case I and II (see text). The thick (green and white) bonds 
show SWBs in H$e- This particular example shows a spin system 
with K = 2, where A* denotes the maximum number of subsystem 
spins that are connected via SWBs with one environment spin (thick 
white lines, see also Section IV). The medium thick (blue) bonds 
show SWBs in H E . 

Afe-1 N E 

H *=- I I I i w>''- do) 

1=1 7=1+1 a=x,y,z 
N S N E 

%=-EI E (id 

(=1 j=l a=x,y,z 

Here, S and / denote the spin-1 /2 operators of the spins of the 
system and the environment, respectively (we use units such 
that h and ks are one). The spin components Sf and /" are 
related to the Pauli spin matrices, for example SJ is a direct 
product of identity matrices and the Pauli spin matrix A a x = 

j ( j q J in position / of the direct product with 1 < i < N$- 

For the geometry of the whole system, we focus on the one- 
dimensional ring consisting of a system with Ns = 4 spin-1 /2 
particles and an environment with Ne spin-1 /2 particles, see 
Fig. Q] Past simulations have shown that a high connectiv- 
ity spin-glass type of environment is extremely efficient to 
decohere a system lfl6l I28l430ll . so we may expect that the 
one-dimensional ring is one of the most difficult geometries 
to obtain decoherence in short times. 

We assume that the spin-spin interaction strengths of the 
system S are isotropic, = / and that only the nearest- 
neighbor interaction strengths £2" ; - and A™, are non-zero. Note 
that for a ring there are only two bonds with strength A", con- 
necting S and E. We distinguish two cases: 

• Case I: The non-zero values of £lf, and A", are gener- 
ated uniformly at random from the range [— £l,£l] and 
[—A, A], respectively. 

• Case II: All non-zero values of the model parameters 
are identical, Clfj = J and A". = /. This corresponds to 
a uniform isotropic Heisenberg model with interaction 
strength /. 



We will see that these two cases show very different scal- 
ing properties of the decoherence depending on the initial 
state. We also investigate the effects of randomly adding small 
world bonds (SWBs) between spins in the system and envi- 
ronment and between spins in the environment (see Fig. [TJ. 

The initial state of the whole system S + E is prepared in 
two different ways, namely: 

• "X": We generate Gaussian random num- 
bers {a(j,p),b(j,p)} and set c(j,p,t = 0) = 

(a(j,p) + ib(j,p))/^j iP (a 2 (j,p) + b 2 (j,p)). 
Clearly this procedure generates a point on the 
hypersphere in the D-dimensional Hilbert space. Alter- 
natively, we generate points in the hypercube by using 
uniform random numbers in the interval [—1,1], Our 
general conclusions do not depend on the procedure 
used (results not shown). 

• UDUDY: The initial state of the whole system is a 
product state of the system and environment. In this 
paper (Ns = 4), we confine the discussion to the state 
UDUDY, which means that the first, second, third, and 
fourth spin are in the up, down, up, and down state re- 
spectively, and the state of the remaining spins is a "Y" 
state in the (D/2 4 ) -dimensional Hilbert space. The "Y" 
state of the environment is prepared in the same way as 
the "X" state of the whole system. 



III. SCALING ANALYSIS OF cr 

All simulations are carried out for a system S consisting 
of four spins (Ns = 4) coupled to an environment E with 
the number of spins Ne ranging from 2 to 30. The interac- 
tion strengths , j with 1 < i < Ns — 1 are always fixed to 
/ = —0.15. For case I all non-zero D."j and A?, are randomly 
generated from the range [—0.2,0.2]. For case II all non-zero 
£2". and A", are equal to / = —0.15 (isotropic Heisenberg 
model). 



A. Verification of scaling: cases I and II with "X" 

We corroborate the scaling property of Eq. ((8} by numeri- 
cally simulating the quantum spin system (see Eq. (O through 
(fTTI)). If we choose the initial state of the whole system to be 
an "X" state, then during the time evolution the whole system 
will remain in the state "X". Hence, the condition to derive 
Eq. ([H) are fulfilled. Fig. [2] demonstrates that the numerical 
results for both cases I and II agree with Eq. dHJ. In par- 
ticular the insets in Fig. [2] show that for both cases I and II, 
ln(2<r) ps —Ne/2, and that a scales as 1/ \fb~E even if Ne = 2 
and N s =4 (Ne <N s ). 
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FIG. 2. Simulation results for o"(?) (see Eq. ®) for case I (top) 
and case II (bottom) for different sizes N = Ne + 4 of the whole 
system. The initial state of the whole system is "A"' (see text). Curves 
from top to bottom correspond to system sizes ranging from N = 6 
to N — 34 in steps of 2. The insets show the time-averaged values 
of (7(f) (pluses) as a function of the size Ne of the environment, 
confirming the theoretical prediction of Eq. l(8) (solid line). 



B. Different initial conditions 

We investigate the effects of the dynamics by preparing the 
initial state of the whole system such that it is slightly differ- 
ent from "X". The initial state of the whole system is set to 
UDUDY. In contrast to Fig. [2] we will see that the two cases I 
and II behave differently. 




1000 2000 3000 4000 5000 6000 7000 8000 
t 



FIG. 3. Simulation results for a(t) (see Eq. lf5J) for case I for dif- 
ferent sizes N = Ne + 4 of the whole system. The initial state of 
the whole system is UDUDY (see text). Curves from top to bottom 
correspond to system sizes ranging from N = 6 to N = 34 in steps 
of 2. The inset shows the time-averaged values of a(t) (pluses) as a 
function of the size Ne of the environment. The data obey the scaling 
property of Eq. (|8) (solid line). 




200 400 600 800 1000 1200 1400 
t 



FIG. 4. Same as Fig.[3]for case II instead of case I. Curves from top 
to bottom correspond to system sizes ranging from N = 16 to N = 34 
in steps of 2. The solid line in the inset is a guide to the eyes. 



2. Case II and UDUDY 



1. Case I and UDUDY 

In Fig. [3] we present the simulation results for case I, the 
couplings in the Hamiltonians He and H$e are chosen uni- 
formly at random. The size N = Ne + 4- of the whole system 
ranges from 6 to 34. An average over the long-time station- 
ary steady-state values of o(t) still obeys the scaling prop- 
erty of Eq. ©, showing that a decreases as 1 / y/Dg, where 
D E = 2 Ne . If N E ->• °°, O ->• 0. This suggests that in the ther- 
modynamical limit the system S decoheres completely. 



We consider the case in which the whole system is de- 
scribed by the isotropic Heisenberg model (J" +1 = = 
A" +1 = J). In Fig. 2] we present simulation results for dif- 
ferent system sizes N = Ne + 4 ranging from 16 to 34. From 
Fig- III it is seen that the behavior for case II is totally differ- 
ent from that of case I (see Fig. |3). In particular, a{f) does 
not scale with the dimension of the environment. From the 
present numerical results, we cannot make any conclusions 
about the limit for large Ne- However if a(t) approaches zero 
as Ne —> °° (see the fifth column of Table [TJ> it does so very 
slowly. 
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FIG. 5. (Color online) Simulation results for cx(f) for case I for N = 
22 and N s = 4. Red solid line: the initial state is UDUDY (see text); 
green dashed line: the initial state is "X" (see text). 
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FIG. 6. (Color online) Difference Act between the time-averaged 
values of a(t) for the initial state UDUDY and "X" of the whole 
system (see Table|I]( as a function of the size of the environment Ng. 
Pluses: case I; circles: case II. The dotted line is a linear fit to the 
data (pluses) for the UDUDY initial state, excluding the first three 
data points, resulting in Ao" = 0.049/ \/De- 



C. Computational effort 



TABLE I. The time average of a(t) in the stationary regime shown 
inFigs.[U[3]and|4] 





prediction 




case I 
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can understand the very different behavior of cases I and II, 
see Figs. [3] and |U by considering the stationary states that are 
obtained. Figure|5]shows that the final values of o(t ) for case I 
are very close for both initial states "X" and UDUDY. This 
suggests that the final stationary state in case I has proper- 
ties similar to those of a state "X", and hence case I obeys 
the scaling property of Eq. (|H) to a good approximation. The 
time-averaged values of a(t) in Figs. [2] and |4] denoted by 
<7, are listed in Table U From Table [[] we see that the values 
of cf for case II with an initial state UDUDY are very dif- 
ferent from those with an initial state "X", and do not show 
the scaling property of Eq. (|H). Thus, the numerical results 
suggest that the initial state and the randomness of the inter- 
action strengths play a very important role in the dynamical 
evolution of the decoherence of a system coupled to an envi- 
ronment. In particular, for case II, starting from a state "X" the 
time-averaged values of <j(f ) scale asc« 1/ y'Dg, but such 
scaling is not observed for starting from a state UDUDY. 



In this paper, the largest number of spins that we simulated 
is N — 34. Using the Chebyshev polynomial algorithm and 
a large time step (t rj IOtt), the N = 34 simulation for the 
bottom curves in Fig.|2](up to a time t w 600) took about 0.3 
million core hours on 16384 BG/P (IBM Blue Gene P) pro- 
cessors, using 1024 GB of memory. Similarly, it took about 4 
million core hours to complete the N = 34 curve in Fig.[3](up 
to a time t w 8000). 



D. Summary: initial state dependence 

For an initial state "X" of the whole system the scaling of <7, 
as given by Eq. ((H), works extremely well for both case I and 
case II, as seen in Fig.|2] When the initial state is UDUDY, we 



From Tabled] it is seen that the values of a for case I with 
the initial state UDUDY are always slightly larger than those 
with the initial state "X". Therefore, it is interesting to exam- 
ine the difference Ac between the values of <7 for the initial 
states UDUDY and "X". Figure |6] shows that Act for case I 
(red pluses) also scales as 1 / \[De (dotted line), except for the 
first three data points, which is probably due to large fluctua- 
tions in the calculations for these small system sizes. There- 
fore, the dynamics of case I will drive the system to a state 
"X" only when the environment approaches infinity. Figure [6] 
also shows that Aa for case II (circles) is almost constant for 
system sizes Af ranging from 16 to 34. Hence, it is unlikely 
that case II with the initial state UDUDY will decohere, even 
if the simulations could be performed for much longer times 
and for larger system sizes. 
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IV. CONNECTIVITY: RING WITH SMALL WORLD 
BONDS 

We investigate the effects of adding small world bonds 
(SWBs) to the Hamiltonians H$e or/and He for both case I and 
case II (see Fig.[T|). To analyze the addition of SWBs to H$e 
we distinguish between spin systems with K < 2 and K > 2, 
where K denotes the maximum number of subsystem spins 
that are connected via SWBs with one environment spin. This 
distinction is motivated by the distinct decoherence character- 
istics for systems with K < 2 and K > 2 for case I (see next 
subsection). An example of a spin configuartion with K = 2 
is shown in Fig.Q] In particular, we are interested in whether 
systems with SWBs will exhibit the same scaling, and whether 
they will decohere from an initial state faster than either of the 
cases studied thus far. The addition of many SWBs changes 
the graph from a one-dimensional ring to a graph with equal 
bond lengths that can only be embedded in high dimensions. 
The initial states are always UDUDY. Furthermore, in order 
not to change too many parameters simultaneously we start all 
simulations from the same state "F" of the environment. Fur- 
thermore, after choosing the random location (and couplings 
Q. a and A" for case I) of the first SWB we preserve this bond 
when adding additional SWBs. We will see that case I and 
case II still behave very differently. 



A. Case I and SWBs 

For investigating the universality of the final value of a{t ) 
we add SWBs (random couplings in the interval [—0.2,0.2]) 
in the Hamiltonian H$e or/and He for case I, and perform sim- 
ulations for N = 24 with N$ = 4. From Fig. |7^, we see that 
adding more and more SWBs to He speeds up the decoher- 
ence process and that the final value of <j(f) corresponds to 
the one given by Eq. ©. As seen in the inset, adding SWBs 
to He has no noticeable effect on the early time behavior of 
ff(t). 

Adding SWBs exclusively to H$e speeds up the deco- 
herence process even further and even at early times clear 
changes in a(t ) can be observed (see Figs. [7J), c). For spin 
configurations with K = 1, c(f) reaches the value given by 
Eq. ([8]) for sufficiently long times, as can be seen from Fig.[7jS. 
However, for configurations with K = 2 (see Fig. [7]:) or K > 2 
(results not shown) a(t) does not obey the scaling property 
Eq. (|8}. Restoring this scaling property seems to require 
an environment that is much more complex than the one- 
dimensional one as indicated by Fig. [7Jl in which we present 
simulation results for the case that SWBs between all non- 
neighboring environment spins have been added. 

B. Case II and SWBs 

For case II, isotropic SWBs are added to H$e or/and He- 
From Fig. [HJ it is clear that even for long times none of the 
curves approach the dotted horizontal line, the value of o(t ) 
for an initial state "X". Adding SWBs exclusively to He does 



not have a dramatic effect on a(t) and has very little effect at 
early times (see Fig. EK). 

Just as for case I, it is seen that adding a few SWBs ex- 
clusively to H$e for a spin configuration with K = 1 signif- 
icantly decreases the time to approach the steady state, and 
that the SWBs in H$e also lead to a decrease in for a 
fixed time even at early times (see Fig. [Hp). For spin con- 
figurations with K = 2 case I and case II seem to have simi- 
lar decoherence properties if SWBs are added exclusively to 
H$e, as seen by comparing Fig.[7J: and Fig.|8};. However, con- 
necting in addition each pair of non-neighboring environment 
spins by isotropic SWBs drives the curves very far away from 
the value of a(t ) for an initial state "X" (see Fig.[HJi). 



C. Summary: SWBs 

Adding SWBs to H$e or/and to He changes the rate of de- 
coherence as seen by the approach to the asymptotic value for 
o(t). In case II, adding isotropic SWBs to H$e or He effec- 
tively alters some spin-spin correlations leading to a decrease 
in the steady-state value of o{t). However, this decrease is 
not sufficient to reach the steady-state value of c(f ) that com- 
plies with the prediction Eq. ([HJ. Adding isotropic SWBs to 
Hse and connecting in addition each pair of non-neighboring 
environment spins by isotropic SWBs drives the curves very 
far away from the value of a(t) for an initial state "X", even 
much further away than the steady-state value for a ring with- 
out SWBs. In contrast to case I systems with K < 2 and K>2 
do not behave significantly different. 

Comparing case II with case I for K < 2, we conclude 
that without introducing the randomness in the x, y, z com- 
ponents of the spin-spin couplings, the dynamics cannot drive 
the system to decoherence if the initial state is different from 
a state "X". Increasing the complexity of the environment by 
adding isotropic SWBs between all non-neighboring environ- 
ment spins does not help in this respect, even on the contrary. 
However, for case I and configurations with K > 2, increasing 
the complexity of the environment by adding SWBs between 
all pairs of non-neighboring environment spins allows the dy- 
namics to drive the system to decoherence. 

For both case I and case II, adding SWBs in Hse and He 
separately speeds up the decoherence in that it evolves more 
quickly to a stationary state. The asymptotic value for a{t) 
is approached much faster when adding SWBs to Hse instead 
of He, and the SWBs in Hse also affect a(t) at early times. 
Thus a random SWB coupling to the system via Hse is the 
most effective way to decrease the time for decoherence. 



V. RANDOMNESS IN THE ENVIRONMENT 

Section llll Al shows that for the initial state "X" the scaling 
predicted by Eq. (|8j is confirmed both for case I and case II 
(see Fig.|2]i. However, section HIl Bl shows that starting from 
the initial state UDUDY this scaling is approached as 1 / \/D~e 
for case I (see Figs. [3]and|6]l but not for case II (see Figs. [4] 
and|6|. SectionllVlshows that adding SWBs in case II does not 
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FIG. 7. (Color online) Simulation results of a(t) for case I with N = 24 and N s = 4 with SWBs added. The initial state is UDUDY. The 
dotted horizontal line represents the value of Eq. ((8}. Red solid line: ring without SWBs. (a) Randomly added SWBs between non-neighboring 
environment spins. Green long-dashed line: one SWB; orange dotted line: two SWBs; purple short-dashed line: four SWBs; blue dotted- 
dashed line: eight SWBs. (b) Randomly added SWBs between the system and environment spins such that K = 1. Green long-dashed line: 
one SWB; orange dotted line: two SWBs; purple short-dashed line: four SWBs; blue dotted-dashed line: eight SWBs. (c) Randomly added 
SWBs between the system and environment spins such that K = 2. Green long-dashed line: two SWBs; orange dotted line: four SWBs; purple 
short-dashed line: six SWBs; blue dotted-dashed line: eight SWBs. (d) Same as (c) except that each pair of non-neighboring environment 
spins is connected by a SWB. Insets: time evolution for short times. 



significantly change the long-time behavior of a(t) approach- 
ing the predicted value of Eq. ((H). Therefore the natural ques- 
tion to ask is how much randomness is required for o(t) to 
obey the scaling relation Eq. ((H). To answer this question, we 
start from the isotropic Heisenberg ring (case II) and replace 
the interaction strengths of a few randomly chosen bonds by 
random £lfj (see Eq. (fTUb). 

Figure [9] presents the simulation results for o(t) by intro- 
ducing 1, 2, 4, 6 and 8 random bonds in the environment 
Hamiltonian He of Eq. ( fTQ) . The interaction strengths £lf- 
of these randomly selected bonds are drawn randomly from 
a uniform distribution in [—0.2,0.2]. Furthermore, the ran- 
domly selected bond for the case with 1 random bond is also 
a random bond for the case with 2 and more randomly chosen 
bonds, thereby not changing too many parameters at a time. 
Simulations up to time t = 6000 show that introducing 4, 6 
and 8 random bonds leads the system to relax to the predicted 
value of a (see Eq. ((H)). For times up to t = 6000 the effect of 
one or two random bonds is not apparent. Therefore for these 
two cases we performed extremely long runs as shown in the 



inset of Fig. |9] The inset shows that even one random bond 
suffices to recover the asymptotic value Eq. ((H). However the 
time scale to reach the asymptotic value of a can become ex- 
tremely long. We leave the question of how fast the approach 
to the predicted value of a is for future study. 

For understanding the behavior of a(t) in case II with ran- 
domness, we investigate the individual components of the re- 
duced density matrix p for the ring system. We study the 
addition of one, two, up to eight randomly replaced bonds in 
the environment. Recall that once the position for one random 
bond is chosen, this is also one of the random bonds when 
there are two or more random bonds. Similarly, the locations 
of the random positions for a large number of random bonds 
include the same positions and strengths as for a smaller num- 
ber of random bonds. Furthermore, the same initial state "Y" 
of the environment is chosen for all simulations. We stud- 
ied the effect of varying the positions of the randomly chosen 
bonds and of different initial states "Y" for the environment 
for a couple systems and did not find significant changes in 
our observations. 
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FIG. 8. (Color online) Simulation results of o"(?) for case II with N = 26 and Ns = 4 with isotropic SWBs added. The initial state is UDUDY. 
The dotted horizontal line represents the value of Eq. (|8}. Red solid line: ring without SWBs. (a) Randomly added SWBs between non- 
neighboring environment spins. Green long-dashed line: two SWBs; orange dotted line: four SWBs; purple short-dashed line: six SWBs; blue 
dotted-dashed line: eight SWBs. (b) Randomly chosen SWBs between the system and environment spins such that K = 1. Green long-dashed 
line: two SWBs; orange dotted line: four SWBs; purple short-dashed line: six SWBs; blue dotted-dashed line: eight SWBs. (c) Randomly 
chosen SWBs between the system and environment spins such that K = 2. Green long-dashed line: two SWBs; orange dotted line: four 
SWBs; purple short-dashed line: six SWBs; blue dotted-dashed line: eight SWBs. (d) Same as (c) except that each pair of non-neighboring 
environment spins is connected by a SWB. Insets: time evolution for short times. 



Figure [10] presents the results of the time evolution of the 
absolute value |p, ; | of the individual components of the re- 
duced density matrix. For completeness we show both the 
diagonal components and the off-diagonal components. Fig- 
ure [TO] shows that most of the 120 off-diagonal components 
quickly relax to a small value (1 14 black lines in Fig. [T0l(b)- 
(e)). The slowest decaying |py| are plotted in red. There 
are six such components. In the steady state all |p~y| oscil- 
late but have nearly the same time-averaged value, in agree- 
ment with the mean-field-type argument given in AppendixlBl 
Thus, only a few |p~y| are responsible for the lack of scaling 
of a in case II when starting from the initial state UDUDY, 
and also for the long times required to approach the predicted 
value of Eq. (O of a{t ) in the case that there are one or two 
random bonds. 



VI. CONCLUSIONS AND DISCUSSION 

The main theoretical result of the current paper is Eq. (|8} for 
the decoherence of a quantum system S coupled to a quantum 
environment E. For studying decoherence we examine a(t ), 
which is the square root of the sum of all the off-diagonal 
elements of the reduced density matrix p for S in the basis 
that diagonalizes the Hamiltonian H$ of the system S. We find 
(see also Eq. dS)) that 

a«-^=(l-J-V (12) 

where the reduced density matrix p for S is a D$ x D$ matrix 
while the density matrix of the whole system S + EisaDxD 
matrix with D = D$De- Thus De does not have to be very 
large in order for the predicted scaling to hold, in particular 
the scaling requires De > 1 > DJ . In addition the scaling 
requires that S + E is driven from an initial wave function to- 
ward a steady state which is well described by a state which 
we called "X". 
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FIG. 9. (Color online) Simulation results of a(t) obtained by se- 
lectively replacing isotropic spin-spin interactions by random bonds. 
The size of the system and whole system are N$ = 4 and N = 26, 
respectively. The initial state is UDUDY. Red solid line: 1 random 
bond; green long-dashed line: 2 random bonds; purple dotted line: 
4 random bonds; orange short-dashed line: 6 random bonds; blue 
dotted-dashed line: 8 random bonds. Inset: simulation results for 
one and two random bonds for long times. 



We have performed large-scale real-time simulations of the 
time-dependent Schrodinger equation for Ns spins in the sys- 
tem and Ne spins in the environment. We have simulated spin- 
1 /2 systems with N = N s + Ne up to N = 34, all with Ns = 4. 
Starting from a state "X" for S + E the simulations agree very 
well with the scaling prediction Eq. ( fT2] >. as shown in Fig. [2] 
In Appendix ICl we demonstrate that in this case not only the 
off-diagonal elements of p obey a scaling relation but also its 
diagonal elements obey a scaling relation, although a different 
one. 

Therefore as long as the dynamics drives the initial state to 
a state "Z" which has similar properties as "X" the scaling re- 
lation Eq. dTzb should hold. The next step is to examine under 
what conditions our test quantum model is driven to the state 
"Z", and study the time scale needed to relax from an initial 
state to the state "Z". For the one-dimensional quantum spin- 
1/2 ring we find that homogeneous couplings do not lead to 
an evolution to the state "Z" (Fig. |4j, and hence the scaling 
as 1 I\[De is not observed. This conclusion is not modified 
if some randomly chosen homogeneous small world bonds 
are added (Fig. [8j. Also systems with random couplings and 
random small world bonds between system and environment 
spins such that the maximum number of system spins that in- 
teract with one environment spin is two or larger do not evolve 
to a state "Z" (Fig.|7J;). In this case, the environment requires 
a more complex connectivity than the simple one-dimensional 
one in order to observe the scaling as 1 / \[De (Fig.[7Jl). There- 
fore, although we find that some randomness in the interaction 
strengths in E or between S and E the dynamics is very im- 
portant to drive the whole system toward the state "Z", as seen 
in Figs.[3][5l|7h>b> and|9]it is not always sufficient. Moreover 
it may take a long time to evolve toward the state "Z" if there 
is only a little randomness (Fig. |9]l or if the environment E is 
large (the N = 3A results of Fig. |3J. The long time that may 
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FIG. 10. (Color online) Time evolution of the components pj; of the 
reduced density matrix of the system with N = 26 and N$ = 4. The 
initial state is UDUDY. (a) Case II. Starting from case II, 1 (b), 2 
(c), 4 (d), 6 (e) and 8 (f) random bonds are introduced in He- Blue 
lines: diagonal components \pa\; red lines: all 6 slowly decaying 
components for |py | for one random bond; black lines: all other 1 14 
off-diagonal components |py|. 



be required to approach the state "Z" is due to only a few off- 
diagonal elements of p, as seen in Fig. |T0j We find that the 
approach to the state "Z" can be sped up by adding random- 
ness to E (Figs.l9landllOI>. 

What do our results say about the approach to the quan- 
tum canonical ensemble? The canonical ensemble is given by 
the diagonal elements of the reduced density matrix p if the 
off-diagonal elements (as measured by a(t )) can be neglected 
01 [l7t] . As long as E has a finite Hilbert space De our scaling 
results can be used to argue that in a strict sense, the system 
will not be in the canonical state unless De — > °°. However, if 
the canonical distribution is to be a good approximation for 
some temperatures T up to some chosen maximum energy 
£hoid > 0, then this requires that exp ( — £hoid/ ksT)^$>a which 
gives for our spin-1 /2 system kgT ^> 2E\ m n/ [A^-ln(2)]. For 
this argument to hold in the canonical distribution the energies 
are taken to be positive values above the ground state energy. 
This lack of thermalization at low temperatures for small sys- 
tems is supported by simulations in Ref . lfl7ll . 

What do our results say about trying to prolong the time 
to decoherence in order to build practical quantum encryption 
or quantum computational devices? The important thing is 
to ensure that the system is not driven toward the state "Z", 
or at least that it takes a very long time to approach the state 
"Z". This can be achieved by changing the Hamiltonian of 
the system, H = H$ + He +H$e, such that it has very small 
randomness particularly in the coupling between the system 
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and the environment, Hse- Alternatively extrapolating from 
Fig. [10] if one can devise an experimental procedure, for ex- 
ample a time-dependent procedure, to keep even a few of the 
off-diagonal elements of p large then the scaling prediction 
Eq. ( TT2b for the decoherence can be avoided, at least for rea- 
sonable timescales. 

The scaling of Eq. ( fT~2l > can be contrasted with the pre- 
dicted scaling of the Hilbert space variant of a whole sys- 
tem which should be proportional to (D + 1 ) ~ 1 for the ex- 
pectation value of a local operator [31]. The results of the 
current research are also relevant for methodologies for me a- 
suring finite-temperature dynamical correlations [1321 without 
performing the complete TDSE evolution of the whole sys- 
tem. 

We leave as future work the coupling between a system S 
composed of spin- 1/2 objects (qubits) and an environment E 
composed of harmonic oscillators. In particular, we have re- 
cently been able to build on exact calculations of a single spin 
coupled to specific types of spin environment ll33tl to devise 



an algorithm that does not have computer memory constraints 
limited by the size of De l34l l35tl . We are working to extend 
this algorithm to other types of environment and for more than 
one spin in the system S. 



Appendix A: Scaling without an environment 

For comparison of the scaling of a{f) for the cases with and 
without an environment, we derive the scaling for the case of 
no environment. In the energy basis |;) of the (system, which 
is now the whole system) Hamiltonian H, the density matrix 
has elements 

p iJ (t)=c l (t)c]{t). (Al) 

We use from Ref. (0 the equations (A. 12) and (A.23). The 
expectation value is 



D s D S \ D S As 



E (2a 2 ) =E £ £ \c t (t)cj(t) \ A \ci(t)cj(t) 



f)|2 M f) = 1 -^ = §Sr (A2) 



The final scaling result for the quantity a that we measure is 



a«_Ly£(2^_-L^_i_-L _J_ + _J_ _L_ + _±_ + .... (A3) 



Therefore without an environment, a approaches a constant 
as the size of the system (which is the whole system) grows. 
This also means that for the state "X", if all off-diagonal 

i 1 2 

elements are the same they will have a size of P(/(f) = 
l/Ds (Ds — 1) ~ l/^l while if all the diagonal elements are 
equal (corresponding to infinite temperature) |p,,(f)| 2 = 1 /Ds 
since Tr p (t ) = 1 . We have performed simulations (results not 
shown) to ensure that for the case without an environment a 
obeys the scaling relation of Eq. ( IA3t and it does. 



Appendix B: Mean-field-like reduced density matrix 

We make a connection between a and the quantum purity 

9* = Tr (^(p) 2 ^j. We assume a 'mean-field-type' structure for 

the reduced density matrix, namely we assume that all off- 
diagonal elements have the same size, e. In our simulations 
we find that in the energy basis the imaginary part of the off- 
diagonal elements are very small, which validates our hypoth- 
esis. However, the signs of the real part of the off-diagonal 
elements are not the same, which brings into question our 
'mean-field-like' assumption. Nevertheless, we make the as- 



sumption that 

/ la 1 
]]D s (Ds-l) 

We introduce the matrix J with all its elements having the 
value 1 , the matrix D which is the diagonal matrix composed 
of the diagonal elements of p, and the identity matrix I. Note 
that J 2 = DsJ- The 'mean-field-type' assumption then reads 

p=D + £j-£l, (B2) 

which as seen from the graphs in Fig. [10] should be a reason- 
able assumption in the steady state regime. We will use the 
relationships 

Tr(D) = l, 
Tr(D 2 ) < 1, 

Tr(I)=Tr(J) =J D 5 , 
Tr(DJ)=Tr(JD) = l, 
Tr(j 2 )=D 2 , (B3) 

with the first relationship being a consequence of the trace of 
a density matrix being equal to unity. Then one has that 

^ = Tr(p 2 ) 
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Tr((D + eJ-£l) 2 ) 



-2cj z 



: Tr (D 2 - 2eD + £ 2 I + £DJ + £JD 2£ 2 J + £ 2 J 2 ) 
:Tr(D 2 ) 

Tr(D 2 ) 

Tr (D 2 ) 

1 f t 1 1 1 

-I 1 1 = 

D E \ D s D E D S D E D\ 



D E + 
1 



Ds 



(B4) 



In the canonical ensemble the diagonal elements of the re- 
duced density matrix are related to the terms in the canonical 
partition function, in particular p„- = e~^ E '/Z fl^[l7ll . There- 
fore we have a connection between the quantum purity and 
how close the system is to a canonical ensemble. In the steady 
state this difference is of the order of 1 /D E . 

With the same 'mean-field-like' assumption for p in the 
steady state one can look at corrections to the von Neumann 
entropy of the system, 5? = — Tr(pTnp). However, we do not 
find the final result too enlightening. 



Appendix C: Diagonal elements of the reduced density matrix 

In the main text, we investigated the scaling property of the 
off-diagonal elements of the reduced density matrix of a sys- 
tem coupled to an environment. For being complete in the 
contents, we present some numerical and analytical results 
concerning the diagonal elements. 

In general, based on the fact that the system decoheres, i.e. 
the off-diagonal elements of the reduced density matrix ap- 
proach zero, we expect that the diagonal elements take (ap- 
proach to) the form of the canonical distribution exp(— fiEj) 
where j3 = 1 /kaT with T denoting the temperature and £g 
Boltzmann's constant, which is taken to be one in this paper, 
and where Efs denote the eigenvalues of Hg Ql [H. The dif- 
ference between the diagonal elements p„- (/) and the canoni- 
cal distribution is conveniently characterized by 



8(t) 



D S ( D S N 



(Cl) 



with a fitting inverse temperature 

L ,./,.. /, Pnpa(0 - \npjj(t)}/(Ej - E,) 



b(t) 



Li<j,EffkEj 1 



(C2) 



If the system relaxes to its canonical distribution both 8(t) and 
a if) are expected to vanish, b(t) converging to the effective 
inverse temperature b. 



The numerical simulations of which we present the results 
correspond to those used to make Fig. [2] The initial state for 
those simulations is "X". We analyze the diagonal elements, 
instead of the off-diagonal elements, of the reduced density 
matrix and calculate the quantity 5(f). In Fig. [TT] we present 
the time-averaged value 8 of 8{t) for each system size. It is 
interesting to see that the quantity 8 also has a kind of scaling 
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FIG. 11. Simulation results for the time-averaged value 8 of 8(t) 
(see Eq. l ICU ) for case I (bullets) and case II (squares) for different 
sizes N of the whole system. The initial state of the whole system is 
"X" (see text). The dotted line is l/VW. 



property. As the whole system size increases, 8 decreases 
as l/\/D, where D = 2 N . 



In fact the fitting inverse temperature b(t) is very close to 
zero for reasonablely large N E (data not shown). The canon- 
ical distribution of S at b = is represented by a diagonal 
density matrix with elements 1/D$, where D$ = 2" s . Then, 
we are able to derive the scaling property for 8 as we did to 
obtain Eq. (0. The expectation value of 8 is given by 



£(<s 2 ) = £|E 



Ds D E 2 2 

H H E [\ C i-p\ \ C i,p'\ 



D E i 
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(C3) 



From Eq. <[C3}, we have 8 « 1 /Vd for D s > 1 and D E > 1. 
Therefore, if the size of the environment goes to infinity with 
the final state being the state "X", the diagonal elements of the 
reduced density matrix of the system approach 1 /D$. 
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